{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "218c5a92",
   "metadata": {},
   "source": [
    "# 习题\n",
    "## 习题15.1\n",
    "![image.png](./images/exercise1.png)\n",
    "### 使用numpy实现奇异值分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b9b367a5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[ 0.45  0.89]\n",
      " [ 0.89 -0.45]]\n",
      "S= [3. 2.]\n",
      "V= [[ 0.75 -0.   -0.67]\n",
      " [ 0.3   0.89  0.33]\n",
      " [ 0.6  -0.45  0.67]]\n",
      "A= [[ 1.  2.  0.]\n",
      " [ 2. -0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[1, 2, 0],\n",
    "              [2, 0, 2]])\n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "\n",
    "# 设置精度为2\n",
    "np.set_printoptions(precision=2, suppress=True)\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)\n",
    "Sigma = np.zeros_like(A, float)\n",
    "# Fill the main diagonal of the given array of any dimensionality.\n",
    "# For an array `a` with ``a.ndim >= 2``, the diagonal is the list of\n",
    "# locations with indices ``a[i, ..., i]`` all identical. This function\n",
    "# modifies the input array in-place, it does not return a value.\n",
    "np.fill_diagonal(Sigma, S)\n",
    "calc = np.dot(np.dot(U, Sigma), V)\n",
    "print(\"A=\", calc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66fcd677",
   "metadata": {},
   "source": [
    "### 自编程实现奇异值分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ddf01a14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[-0.45 -0.89]\n",
      " [-0.89  0.45]]\n",
      "S= [[3. 0. 0.]\n",
      " [0. 2. 0.]]\n",
      "V= [[-0.75 -0.   -0.67]\n",
      " [-0.3  -0.89  0.33]\n",
      " [-0.6   0.45  0.67]]\n",
      "A= [[ 1.  2. -0.]\n",
      " [ 2. -0.  2.]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import null_space\n",
    "\n",
    "\n",
    "def my_svd(A):\n",
    "    m = A.shape[0]\n",
    "\n",
    "    # (1) 计算对称矩阵 A^T A 的特征值与特征向量，\n",
    "    W = np.dot(A.T, A)\n",
    "    # 返回的特征值lambda_value是升序的，特征向量V是单位化的特征向量\n",
    "    lambda_value, V = np.linalg.eigh(W)\n",
    "    # 并按特征值从大到小排列\n",
    "    lambda_value = lambda_value[::-1]\n",
    "    lambda_value = lambda_value[lambda_value > 0]\n",
    "    # (2)计算n阶正价矩阵V\n",
    "    V = V[:, -1::-1]\n",
    "\n",
    "    # (3) 求 m * n 对角矩阵\n",
    "    sigma = np.sqrt(lambda_value)\n",
    "    \n",
    "    # np.diag(sigma)：这个函数调用创建了一个对角矩阵，其中 sigma 是一个数组或列表，其元素被放置在对角线上。\n",
    "    # 如果 sigma 是一个标量值，则结果是对角线上所有元素都为该值的对角矩阵。\n",
    "    # np.eye(*A.shape)：这个函数调用创建了一个与矩阵 A 同维度的单位矩阵。\n",
    "    # *A.shape 是将 A 的形状（行数和列数）作为参数传递给 np.eye 函数，这样创建的单位矩阵的维度和 A 是相同的。\n",
    "    # @：这是 Python 中的矩阵乘法运算符，用于执行两个矩阵的乘法。\n",
    "    S = np.diag(sigma) @ np.eye(*A.shape)\n",
    "\n",
    "    # (4.1) 求A的前r个正奇异值\n",
    "    r = np.linalg.matrix_rank(A)\n",
    "    U1 = np.hstack([(np.dot(A, V[:, i]) / sigma[i])[:, np.newaxis] for i in range(r)])\n",
    "    # (4.2) 求A^T的零空间的一组标准正交基\n",
    "    U = U1\n",
    "    if r < m:\n",
    "        U2 = null_space(A.T)\n",
    "        U2 = U2[:, r:]\n",
    "        U = np.hstack([U, U2])\n",
    "\n",
    "    return U, S, V\n",
    "\n",
    "\n",
    "A = np.array([[1, 2, 0],\n",
    "              [2, 0, 2]])\n",
    "\n",
    "np.set_printoptions(precision=2, suppress=True)\n",
    "\n",
    "U, S, V = my_svd(A)\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V)\n",
    "calc = np.dot(np.dot(U, S), V.T)\n",
    "print(\"A=\", calc)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6e8dc3d7",
   "metadata": {},
   "source": [
    "## 习题15.2\n",
    "![image.png](./images/exercise2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4985f4b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[-0.82 -0.58  0.    0.  ]\n",
      " [-0.58  0.82  0.    0.  ]\n",
      " [ 0.    0.    1.    0.  ]\n",
      " [ 0.    0.    0.    1.  ]]\n",
      "S= [5.46 0.37]\n",
      "V= [[-0.4  -0.91]\n",
      " [-0.91  0.4 ]]\n",
      "A= [[ 2.  4.]\n",
      " [ 1.  3.]\n",
      " [-0.  0.]\n",
      " [-0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "# 计算矩阵𝐴A的奇异值分解\n",
    "A = np.array([[2, 4],\n",
    "              [1, 3],\n",
    "              [0, 0],\n",
    "              [0, 0]])\n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "np.set_printoptions()\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)\n",
    "\n",
    "# 根据奇异值分解的结果，写出外积展开式\n",
    "calc = S[0] * np.outer(U[:, 0], V[:, 0]) + S[1] * np.outer(U[:, 1], V[:, 1])\n",
    "print(\"A=\", calc)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c65789b9",
   "metadata": {},
   "source": [
    "## 习题15.5\n",
    "![image.png](./images/exercise5.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1a6e4952",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U= [[ 1.   -0.    0.   -0.01]\n",
      " [ 0.    1.    0.   -0.  ]\n",
      " [ 0.    0.    1.    0.  ]\n",
      " [ 0.01  0.    0.    1.  ]]\n",
      "S= [20.62 10.44  1.    0.97]\n",
      "V= [[ 0.    0.96 -0.   -0.    0.29]\n",
      " [ 0.97 -0.   -0.   -0.24 -0.  ]\n",
      " [ 0.24  0.    0.    0.97  0.  ]\n",
      " [ 0.    0.29  0.    0.   -0.96]\n",
      " [ 0.    0.    1.    0.    0.  ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[0, 20, 5, 0, 0],\n",
    "              [10, 0, 0, 3, 0],\n",
    "              [0, 0, 0, 0, 1],\n",
    "              [0, 0, 1, 0, 0]]) \n",
    "\n",
    "# 调用numpy的svd方法\n",
    "U, S, V = np.linalg.svd(A)\n",
    "\n",
    "print(\"U=\", U)\n",
    "print(\"S=\", S)\n",
    "print(\"V=\", V.T)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2da70967",
   "metadata": {},
   "source": [
    "![image.png](./images/exercise5_2.png)\n",
    "# 自编程实现奇异值分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "16e8f9db",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_U = \n",
      " [[ 0.2  -0.51 -0.  ]\n",
      " [ 0.03 -0.8   0.  ]\n",
      " [ 0.01 -0.27  0.  ]\n",
      " [ 0.18  0.02 -0.  ]\n",
      " [ 0.35  0.05 -0.  ]\n",
      " [ 0.88  0.12 -0.  ]\n",
      " [ 0.18  0.02 -0.  ]]\n",
      "sigmasMatrix = \n",
      " [[9.82 0.   0.  ]\n",
      " [0.   5.26 0.  ]\n",
      " [0.   0.   0.  ]]\n",
      "X_V = \n",
      " [[ 0.58  0.04  0.82]\n",
      " [ 0.58  0.04 -0.41]\n",
      " [ 0.58  0.04 -0.41]\n",
      " [ 0.05 -0.71  0.  ]\n",
      " [ 0.05 -0.71  0.  ]]\n"
     ]
    }
   ],
   "source": [
    "#基于矩阵分解的结果，复原矩阵\n",
    "def rebuildMatrix(U, sigma, V):\n",
    "    a = np.dot(U, sigma)\n",
    "    a = np.dot(a, np.transpose(V))\n",
    "    return a\n",
    "\n",
    "\n",
    "#基于特征值的大小，对特征值以及特征向量进行排序。倒序排列\n",
    "def sortByEigenValue(Eigenvalues, EigenVectors):\n",
    "    # 这是 NumPy 中的一个函数，用于获取数组中元素的排序索引。np.argsort() 返回的是一个索引数组，这个索引数组表示原始数组元素需要如何重新排列才能得到一个排序后的数组。\n",
    "    index = np.argsort(-1 * Eigenvalues)\n",
    "    Eigenvalues = Eigenvalues[index]\n",
    "    EigenVectors = EigenVectors[:, index]\n",
    "    return Eigenvalues, EigenVectors\n",
    "\n",
    "\n",
    "#对一个矩阵进行奇异值分解\n",
    "def SVD(matrixA, NumOfLeft=None):\n",
    "    #NumOfLeft是要保留的奇异值的个数，也就是中间那个方阵的宽度\n",
    "    #首先求transpose(A)*A\n",
    "    matrixAT_matrixA = np.dot(np.transpose(matrixA), matrixA)\n",
    "    #然后求右奇异向量\n",
    "    lambda_V, X_V = np.linalg.eig(matrixAT_matrixA)\n",
    "    lambda_V, X_V = sortByEigenValue(lambda_V, X_V)\n",
    "    #求奇异值\n",
    "    sigmas = lambda_V\n",
    "    sigmas = list(map(lambda x: np.sqrt(x)\n",
    "                      if x > 0 else 0, sigmas))  #python里很小的数有时候是负数\n",
    "    sigmas = np.array(sigmas)\n",
    "    sigmasMatrix = np.diag(sigmas)\n",
    "    if NumOfLeft == None:\n",
    "        rankOfSigmasMatrix = len(list(filter(lambda x: x > 0,\n",
    "                                             sigmas)))  #大于0的特征值的个数\n",
    "    else:\n",
    "        rankOfSigmasMatrix = NumOfLeft\n",
    "    sigmasMatrix = sigmasMatrix[0:rankOfSigmasMatrix, :]  #特征值为0的奇异值就不要了\n",
    "\n",
    "    #计算右奇异向量\n",
    "    X_U = np.zeros(\n",
    "        (matrixA.shape[0], rankOfSigmasMatrix))  #初始化一个右奇异向量矩阵，这里直接进行裁剪\n",
    "    for i in range(rankOfSigmasMatrix):\n",
    "        X_U[:, i] = np.transpose(np.dot(matrixA, X_V[:, i]) / sigmas[i])\n",
    "\n",
    "    #对右奇异向量和奇异值矩阵进行裁剪\n",
    "    X_V = X_V[:, 0:NumOfLeft]\n",
    "    sigmasMatrix = sigmasMatrix[0:rankOfSigmasMatrix, 0:rankOfSigmasMatrix]\n",
    "    #print(rebuildMatrix(X_U, sigmasMatrix, X_V))\n",
    "\n",
    "    return X_U, sigmasMatrix, X_V\n",
    "\n",
    "A = np.array([[1, 1, 1, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0],\n",
    "              [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]])\n",
    "X_U, sigmasMatrix, X_V = SVD(A, NumOfLeft=3)\n",
    "print('X_U = \\n', X_U)\n",
    "print('sigmasMatrix = \\n', sigmasMatrix)\n",
    "print('X_V = \\n', X_V)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
